library(vegan)

# input for vegan

veg_birds <- abn%>%
  filter(Taxa == "Birds")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

veg_butterfly <- abn%>%
  filter(Taxa == "Butterflies")%>%
  group_by(Study, Transect, Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  mutate(sample = paste(Study, "_", Transect))%>%
  group_by(sample)%>%
  #mutate(N = sum(n, na.rm = T))%>%
  #filter(N > 0)%>%#removing samples where no species were detected
  dplyr::select(sample, Specific.Epithet, n)%>%
  ungroup()%>%
  spread(Specific.Epithet, n, fill = c(n = 0))%>%
  column_to_rownames(var = "sample")

# sample info

samp <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  dplyr::select(Study, Transect, Habitat)%>%
  distinct()%>%
  mutate(sample = paste(Study, "_", Transect))

## Species accumulation curves

rarebutts <- rarecurve(veg_butterfly) ## butterflies

# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html

names(rarebutts) <- samp$sample

protox <- mapply(
  FUN = function(x, y) {
    mydf <- as.data.frame(x)
    colnames(mydf) <- "value"
    mydf$sites <- y
    mydf$subsample <- attr(x, "Subsample")
    mydf}, 
  x = rarebutts, y = as.list(names(rarebutts)), SIMPLIFY = FALSE)

rarebutts.df <- do.call(rbind, protox)

rownames(rarebutts.df) <- NULL  # pretty

rarebutts.df$Taxa <- "Butterflies"

## repeating for birds

rarebirds <- rarecurve(veg_birds)

# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html

names(rarebirds) <- samp$sample

protox <- mapply(
  FUN = function(x, y) {
    mydf <- as.data.frame(x)
    colnames(mydf) <- "value"
    mydf$sites <- y
    mydf$subsample <- attr(x, "Subsample")
    mydf}, 
  x = rarebirds, y = as.list(names(rarebirds)), SIMPLIFY = FALSE)

rarebirds.df <- do.call(rbind, protox)

rownames(rarebirds.df) <- NULL  # pretty

rarebirds.df$Taxa <- "Birds"

## joining rarefaction values

rarefied <- full_join(rarebirds.df, rarebutts.df)

## adding sample info

rarefied <- left_join(rarefied, samp, by = c("sites" = "sample"))

# Plot.
ggplot(rarefied, aes(x = subsample, y = value, col = Study))+
  geom_line(size = 1)+
  labs(x = "Sample Size", y = "Species")+
  facet_grid(Taxa ~ Transect)

ggsave(last_plot(), filename = "./Figures and Tables/figureA.png", height = 8, width = 16, dpi = 300)  

# Function to caluculate Shannon diversity

shanon <- function(c){
  
  if(sum(c) < 20){
    
    H <- "Sample size too small"
    
  }else{
    
    require(SpadeR)
    require(tidyr)
    
    div <- Diversity(c, datatype = "abundance")
    
    H <- data.frame(div$Shannon_diversity)%>%
      rownames_to_column("Method")
    
  }
  
  return(H)
  
}

# Calculating diversity in each site over years

div_year <- abn%>%
  group_by(Study, Taxa, Transect,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect)%>%
  nest()%>%
  mutate(H = map(data, ~shanon(.)))%>%
  dplyr::select(-data)%>%
  unnest(H)%>%
  filter(grepl("MLE", Method) | is.na(Method))

# Summary

div_year%>%
  group_by(Taxa, Study)%>%
  skimr::skim(Estimate)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

# t test

div_year%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(Estimate ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()

# plot

div_year%>%
  ggplot(aes(Taxa, Estimate, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Effective number of species")+
  scale_fill_brewer(palette = "Greys")

ggsave(last_plot(), filename = "./Figures and Tables/figure2.png", height = 6, width = 8, dpi = 300)  

# compiling trophic guild data

host <- read.csv("./Data/Butterflies_host plant_280820.csv") #butterfly host plant data

diet <- read.csv("./Data/Birds_diet.csv") #bird diet data

trophic <- host%>%
  full_join(diet)%>%
  dplyr::select(Specific.Epithet, Guild)%>%
  distinct()

# Calculating diversity of guilds

guilds <- abn%>%
  left_join(trophic, by = c("Specific.Epithet"))%>%
  filter(!Guild == "")%>%
  group_by(Study, Taxa, Transect, Guild,  Specific.Epithet)%>%
  summarise(n = sum(Abundance, na.rm = T))%>%
  spread(Specific.Epithet, n, fill = c(0))%>%
  group_by(Study, Taxa, Transect, Guild)%>%
  nest()%>%
  mutate(H = map(data, ~shanon(.)))%>%
  dplyr::select(-data)%>%
  unnest(H)%>%
  filter(grepl("MLE", Method) | is.na(Method))

# summary

guilds%>%
  group_by(Taxa, Guild, Study)%>%
  skimr::skim(Estimate)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Guild, Study, mean, sd)

# t - test

guilds%>%
  filter(Guild != "Liana Specialist")%>% #removing due to lack of obs
  complete(nesting(Study, Taxa), Transect, Guild, fill = list(H = 0))%>%
  group_by(Taxa, Guild)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(Estimate ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()

# plot
guilds%>%
  ggplot(aes(reorder(Guild, Estimate), Estimate, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5, position = position_dodge(preserve = "single"))+
  labs(y = "Effective number of species", x = "Trophic Guild")+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  scale_fill_brewer(palette = "Greys")+
  facet_wrap(~Taxa, scales = "free_x", ncol = 1)

ggsave(last_plot(), filename = "./Figures and Tables/figure4.png", height = 9, width = 8, dpi = 300)  

# trophic specialisation based on Juliard et al. 2006

butterflies_TSI <- host%>%
  mutate(H = length(unique(Hostplant.Family)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Hostplant.Family)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

birds_TSI <- diet%>%
  mutate(H = length(unique(Prey)))%>%
  group_by(Specific.Epithet)%>%
  mutate(h = length(unique(Prey)))%>%
  ungroup()%>%
  mutate(TSI = sqrt((H/h)-1))

CSI.T <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(butterflies_TSI, by = c("Specific.Epithet"))%>%
  left_join(birds_TSI, by = c("Specific.Epithet"))%>%
  mutate(TSI = ifelse(is.na(TSI.x), TSI.y, TSI.x))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, TSI)%>%
  summarise(CSI.T = mean(TSI, na.rm = T))

# summary
CSI.T%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.T)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

# t - test

CSI.T%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.T ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()

CSI.T%>%
  ggplot(aes(Taxa, CSI.T, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Trophic)")+
  scale_fill_brewer(palette = "Greys")+
  scale_y_sqrt()

HSI <- abn%>%
  left_join(transects, by = c("Transect"))%>%
  group_by(Study, Taxa, Habitat, Specific.Epithet)%>%
  # calculating abundance of each species in each habitat type
  summarise(n = sum(Abundance, na.rm = T))%>%
  group_by(Study, Taxa)%>%
  mutate(N = sum(n), # total number of indviduals encountered 
         rel.prop = n/N)%>% # relative proportion of species
  group_by(Study, Specific.Epithet)%>%
  summarise(HSI = sd(rel.prop)/mean(rel.prop))%>%
  replace_na(replace = list(HSI = 0))

CSI.H <- abn%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
  group_by(Study, Taxa, Transect)%>%
  distinct()%>%
  left_join(HSI, by = c("Study", "Specific.Epithet"))%>%
  dplyr::select(Study, Transect, Taxa, Specific.Epithet, HSI)%>%
  summarise(CSI.H = mean(HSI, na.rm = T))

# summary
CSI.H%>%
  group_by(Taxa, Study)%>%
  skimr::skim(CSI.H)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Taxa, Study, mean, sd)

# t - test

CSI.H%>%
  group_by(Taxa)%>%
  nest()%>%
  mutate(test = map(data, ~t.test(CSI.H ~ Study, data = .)),
         sumry = map(test, broom::tidy))%>%
  dplyr::select(sumry)%>%
  unnest()%>%
  tstats()

CSI.H%>%
  ggplot(aes(Taxa, CSI.H, fill = Study))+
  geom_boxplot(width = 0.25, alpha = 0.5)+
  labs(y = "Community specialisation index \n (Habitat)")+
  scale_fill_brewer(palette = "Greys")